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OO , We present here a formulation of the electronic ground-state energy in terms of the second order 



reduced density matrix, using a duality argument. It is shown that the computation of the ground- 
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• the dual cone of A^-representability conditions. Some numerical results validate the approach, both 
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, for equilibrium geometries and for the dissociation curve of N2. 
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I. INTRODUCTION 



As early as in 1951, it was noticed by Coleman that the electronic A^-body ground-state 
energy could be obtained by minimizing over the set of A^-representable two-body reduced 
density matrices (2-RDM), and Mayer definitely opened the field in 1955 with his pioneer- 
ing article At a conference in 1959, Coulson then proposed to completely eliminate 
wavefunctions from Quantum Chemisty, since all the electronic ground-state properties of 
molecular systems can be computed from the 2-RDM Unfortunately, the set of 

A^-representable 2-RDM is not known explicitly. Some mathematical characterizations were 
provided 0,|5, ^ but they could not be used to derive a numerical method with a complexity 
of a lower order than the usual A^-body problem. Analytical approaches for model systems 
(see e.g. 0) were also proposed in order to precise the accuracy of the A^-representability 
conditions in specific cases. In practice, only approximate RDM minimization problems, in 
which only a few necessary A^-representability conditions are imposed (see the geometric 
constraints of 0], or the so-called P,Q,G conditions 0, [3]), can be considered. The first 
numerical studies relying on this strategy gave encouraging results lll |. 

Recently a new interest in the Reduced Density Matrix (RDM) approach arose. Impres- 
sive numerical results have been obtained by two different strategies issued from semidef- 
inite programming: primal-dual interior point methods jl^ . 13, Q, on the one hand, 
augmented Lagrangian formulations using matrix factorizations of the 2-RDM 16. IitI lisj] 
on the other hand. These results use a small number of known necessary conditions of 
A^-representability. Yet, the so-obtained ground-state energies are as accurate as the ones 
obtained with coupled-cluster methods, see e.g. Q^Il^- addition, these energies provide 
lower bounds of the Full CI energies, whereas the variational post Hartree-Fock methods, 
such as CI or MCSCF, all provide upper bounds. 

Although the current implementations of variational 2-RDM algorithms are limited to 
the simulation of small molecules in small basis sets, we believe that improvements of the 
algorithms and increase of computational power will make it possible to simulate larger 
molecules and to use larger basis sets in a near future. This will allow in particular to assess 
the convergence of the RDM approach with respect to the size of the basis set, for a given 
molecular system. 

Since the RDM method is a linear minimization problem over a convex set of complicated 
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structure, it is natural to use the concept of duality to mathematically characterize and 
numerically compute the minimum. Duality is an underlying issue in all the RDM studies 



mm 

prob 



19l . |20| , but surprisingly, the specific form of the dual formulation of the RDM 



e.g 



em has not yet been used to derive an efficient algorithm. The current methods (see, 

nn 
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all use general duality considerations in their algorithms, but none 
of them solves directly (and only) the dual RDM problem. The purpose of the present article 
is to present such an approach. As will be shown below, the associated dual optimization 
problem boils down to the search of the zero of a one-dimensional convex function. 

The paper is organized as follows. After setting the problem in section m we derive 
the RDM and the approximate RDM dual problems by standard Lagrangian methods in 
section lllll Then, in section IIVI we propose a new algorithm which aims at solving directly 
the dual problem. Section |3 eventually presents some numerical results demonstrating that 
this new method is an interesting and efficient alternative to the existing methods. 

II. NOTATION 

Let us consider a finite-dimensional space f) := span(xj, i = l,...,r) where {xi)i>i is a 
Hilbert basis of the one-body space L^(]R^ x{|'f),||)},C). Most of our analysis is also valid 
in infinite dimension but for the sake of simplicity, we restrict to the finite-dimensional case. 
The electronic Hamiltonian acts on the A^-body fermionic space A!)Li ^ antisymmetric 
A^-body wavefunctions "^{xi, ...,xn) and is formally defined as 

^ 1 

1=1 l<i<j<N ' * 

where h = —A/2+V and V is the external Coulomb potential generated by the nuclei. In the 
whole paper, we denote by x = (x, a) the vector containing both the space variable x G 
and the spin variable a G {| T); I i)}- For any vector space X, we denote by S{X) the space 
of self-adjoint matrices acting on X, and by V{X) C S{X) the cone of positive semi-definite 
matrices. We also use the simplified notation Vn ■= V (^Af^ "^n ■= S ^Af^ The 

ground-state energy then reads 

E= inf (*,i/^$)= inf tr(i/jvT). (1) 

1^1=1 tr(T)=l 
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The second equality holds true for the infimum of the energy over the set of mixed states 
coincides with the infimum of the energy over the set of pure states. In mathematical words, 
the minimum of a linear function over a convex set is attained on an extremal point of the 
convex set. 

The 2-RDM F associated with an A^-body density matrix T G Vn is defined by means of 
Kummer's contraction operator as jj, 0] 

n^;^ = ^^(T)^;;^ = iv(iv - 1) ^ T^S:.£. (2) 

fc3,...,fejv = l 

Then, the cone Cn of N-representable two-body density matrices is by definition the image 
by Lj^ of the cone Vn of A^-body density matrices: 

Of course the 2-RDMs of physical interest are the elements F G Cat which arise from a 
normalized A^-body density matrix T, i.e. which additionally satisfy tr(F) = A^(A^ — 1). 

Since the Hamiltonian only contains two-body interactions, the energy, of the system 
can be expressed in terms of the two-body density matrix F only (see, e.g. 



E = inf tT(KNT) 

reCiv, 
tT{r)=N{N-i) 



(3) 



where we have introduced 



hxi + 



2(iV-l) 2|xi-X2r 
Formula © is an obvious consequence of the identity Hn = {Lj^)*KN where (Lj^)* is the 
adjoint of sometimes also called a lifting operator. Notice that we did not impose any 
constraint on the spin state in (jH)), but such constraints can be easily taken into account. 



III. DUAL FORMULATION OF THE RDM MINIMIZATION PROBLEM 

We now present the dual formulation of the minimization We recall that the polar 
cone C* of a cone C in any Hermitian space is defined as C* = {x | Vy G C, {x, y) > 0}, where 
(■, ■) denotes the considered scalar product. The dual method then consists in formulating 
in terms of (Cat)* instead of Cn: 



E = N{N - 1) sup{/i I i^^ - ^ G {Cn)*}. 



(4) 
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We therefore obtain an optimization problem in dimension 1 over /i G M which is the 
variable dual to the constraint tr(r) = A^(A^ ^1)- Of course characterizing the polar cone 
(Ctv)* is as difficult as characterizing C^, this issue is called the N -representability problem. 
Indeed = (Cat)**. Even if the dual formulation (0)) does not simplify the theoretical 
iV-representability problem, it turns out to be more convenient for numerical purposes, as 
will be shown below. 

Formula (j3)) can be easily derived from Q. Introducing the Lagrangian 

£(r, B, /i) = tr(i^7vr) - tr(5r) - /i{tr(r) - N{N - 1)}, 

it follows 

E= inf sup C{V,B,jj). (5) 

It then suffices to exchange the inf and the sup in (0) to obtain (@]). Indeed, it is a general 
fact that for any cone C in a finite-dimensional space 

inf (a, x) = supjyU | a — bfi E C*}. (6) 

xGC, (b,x)=l 



Note that this property has been already used in the RDM setting by Erdahl j20| . We shall 
use it again below. 

Since both (Cn)* and Cn are unknown and difficult to characterize, it is necessary to 
approximate (0)) by a variational problem that can be carried out numerically. To this end, 
some necessary conditions for A^-represent ability are selected. We consider in this paper L 
conditions of the following general form 

V£=1...L, £^(r)>0 (7) 

where for any i, Ce : S2 ^ S{Xi) is a linear map and Xi is some vector space. For instance, 
the so-called P-condition Ci{r) = F > (with Xi = i) A i)) originates from the Kummer 
operator preserving positivity, and will always be considered. Other classical necessary 
conditions of A^-representability will be introduced below. Imposing only the necessary 
conditions ((Tj) means that Cn is replaced by the approximate cone Capp D defined as 

Capp := {TeS2\W = £,(F) > 0}. 

Its polar cone can easily be shown to be 

(Capp)* := I Be e S{Xe), 5, > oj , (8) 



and the associated approximate energy is then, in view of ©, 

Eapp = ^inf tr(K^r) (9) 

tT:{T)=N{N-l) 

= N{N - 1) SUp{/i \KN-f^e (Capp)*}. (10) 

Let us emphasize that, since Capp D C^, the energy i?app is a lower bound to the full 
CI energy in the chosen basis, -Eapp < E. We present below an algorithm for solving 
problem (fTUj) . Notice that we obtain only the ground-state energy (and not the ground state 
density matrix), but, resorting to first order perturbation theory, any observable including 
at most two-body interaction terms can be obtained by a finite difference of energies. 

Some well-known necessary conditions of the form ((7j) are the P, Q, G conditions P,! 
Additional necessary conditions can be considered, such as Erdahl's Ti and T2 conditions 
The P, Q and G conditions correspond to the following linear operators in ((7j): 



Ci{T) = r, 
[^c(r)]-- = -ri- + 5-7^!^, 

where 7^ = -^—^ X]fc=i t is the one-body RDM associated with the two-body RDM F. 
Expressions for the adjoint operators (>Cq)* and (/^c)* were presented in for example. 
Notice that for any F G 1S2, '^^p(F) and >Cq(F) also are antisymmetric, whereas >Cg(F) is not. 
Therefore, Xp = Xq = f) A f) and Xq = t) ® t) in the above general formalism. Notice also 
that Erdahl's three-index conditions Ti, T2 require X^^ = f) ® I) Cg) f). 

Our numerical tests were performed using the P, Q, G conditions but our algorithm for 
solving (fTUj) is valid for any set of necessary conditions of the form (|7j). 

IV. ALGORITHM FOR SOLVING THE DUAL PROBLEM 

Let us introduce the distance to the dual cone (Capp)* 

S{ii) = dist (Kn - /i, (Capp)*) . 
Denoting /i*pp = E^^pp/ {N{N — 1)), the function 6 satisfies the following properties: 
(z) 6 = on (— oo,/i*pp] and is increasing on [/i*pp, 00); 
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(ii) S is convex on R; 

{Hi) 6^ is continuously different iable on M, thus 6 is continuously differentiable on M\ {/i*pp} 
and 



v/i>/i:pp, s'{f,) = - " 7 (11) 

where denotes the projection of K]\r — ji onto the polar cone (Capp)*. 



2l|. To prove (i), one notices that when /i < /i* 



Proofs for [ii) — (in) can be found 
Kn — = -f^AT — /^* + (/^* — /^) belongs to (Capp)* since fi* — E V2 C (Capp)*. To illustrate 
the above properties, we provide a plot of 6{fi) for N2 in a ST0-6G basis set, see Figure ^ 

Figure I 

In order to compute /U*pp, we use a Newton-like scheme that strongly exploits the above 
mentioned properties in a natural way: starting from an initial energy above /i*pp (such as 
the Hartree-Fock energy for instance) and using the convexity of the function 6, the Newton 
algorithm ensures that the energy fi decreases at each step of the optimization process and 
converges to /iapp- The right derivative of 6 at /i*pp being always positive, the convergence 
rate is guaranteed to be at least superlinear. 

Of course, the most difficult part of the algorithm is the computation of the distance 5(/x) 
to the cone, and of the projection of Kn — fi. To this end, we chose to minimize, for a 
given /i, the objective function 

L 2 
/f^-/i-^(/:,)*5, 



MB) = \ 



under the constraints i?^ > (£ = 1...L), according to the definition (jH|) of the polar 
cone (Capp)*. The above minimization is performed using a classical limited-memory BFGS 
algorithm [2^, keeping the last m = 3 descent directions. The positivity constraints were 



16 



3- 



parametrized by = (C^)^ with Ce symmetric, as suggested by Mazziotti in 

Computing 6{fi) with sufficient accuracy when fi is close to yU*pp can be difficult because 
the minimization of J^{B) then is ill-conditioned. We therefore consider a "truncated" 
version of the Newton algorithm where is updated by a fraction < a < 1 of the Newton 
step. We then use the linearity of 5 for values close to /i*pp to devise a stopping criterion 
limiting the number of iterations. The algorithm is as follows: 
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Algorithm 1 Consider an initial value fi'^ (for example the Hartree-Fock value //hf ), o-nd 
< a < 1. Compute the projection A^o of Kn — /x*^ on (Capp)* and the distance = S{fi^), 
and consider /x^ = /i° — yg-^- For n > 1, 

• Step 1. Compute the projection = ^^^i{Ci)* [(C")^] of — /i" on (Capp)*, the 
associated distance rf" = = \ \Kn — /i" — A^u\\ and the derivative 

• Step 2. Compute the interpolation slope p^ 



• Step 3. If p"' < (1 + e)5'(yu"), then the linear assumption is satisfied and the final value 
is extrapolated from the current position as /i* = — jrj-^; 

• Step 4- Otherwise, set = /i" — ay^-^ and start again from (1) using as initial 
guess = €"1 for any 1=1.. .L. 

In practice, tlie above algoritiim converges in a few iterations. The only time consuming 
step is the projection performed in Step 1. As described above, this projection is done 
iteratively by minimizing the objective function by a limited-memory BFGS algorithm. 
The cost of one BFGS iteration scales as O(r^). We did not observe a clear scaling of the 
number of BFGS iterations with respect to the basis set size. The memory requirements 
scale as O(r^). Both computational time and memory requirements are comparable to those 

of H- 



V. NUMERICAL RESULTS 



We have tested the method on several molecules at equilibrium geometries using data 



from 



for ST0-6G and 6-3 IG basis sets. The results are reported in Table H] and HTl 



respectively. 

Tables I and II 



The reference Full CI (FCI) energies have been computed using GAMESS 2J|. The 



correlation energies are recovered with a good accuracy This is consistent with previous 
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15 



16|,|l7|. 



results already obtained with different RDM methods 

In general, we have observed that the function 5 is almost linear in quite large a right 
neighborhood of /i*pp (see FigureHJ. One iteration of the Newton algorithm already provides 
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a very correct approximation of the exact RDM energy, even when starting from the Hartree- 
Fock level. Usually, only 3 or 4 Newton iterations are necessary to achieve convergence. 
Therefore, the only limiting step of the method is the computation of the distance S{fi) and 
of the projection of K^^ — /i on the polar cone. The method is very robust with respect 
to initial choices of the energy fi^ and the matrices C^. However, we have observed that 
the computational time needed for finding the projection highly depends on the quality 
of the initial guess. The choice of genuine initial conditions is not obvious since we are 
manipulating abstract objects (dual elements of 2-RDM). Some CPU times are reported in 
Table ITTTl for very crude initial conditions = Id and fj,^ ^ O.Q/^hf- 

Table III 

We would like to underline that our projection algorithm is far from being optimal. There 
is clearly much room for improvement here. Let us also mention that the curve fi — 
can be easily sampled using parallel computing (one value of // per processor). 

We also present in Figure |21 dissociation curves for N2 in a ST0-6G basis set. This 



example was already studied in several works j25 



26 



27|. The agreement of our results with 



the reference Full CI is excellent, and the dissociation energy is therefore recovered with a 
very good accuracy. 

Figure II 
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Table and Figures captions 

• Figure I. Left: Distance S{fi) of Kn — /i to the cone (Capp)* as a function of fi for N2 
in a ST0-6G basis set. The tangent at the estimated value for yU*pp is also displayed 
(dotted line). Right: Zoom near the FCI reference value. The Hartree-Fock value is 
/iHF = -1.4435153 while the reference FCI value is /^ci = -1.4453909. 

• Table I. Correlation energies in a ST0-6G basis set. 

• Table II. Correlation energies in a 6-3 IG basis set. 

• Table III. CPU time (s) in a ST0-6G basis using very crude initial guesses {Ci = I). 

• Figure II. Dissociation curve for N2 in a ST0-6G basis set. 
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TABLE I: Cances et al., Journal of Chemical Physics. 



Oy o Lclll 


r ^1 ciici g,y 


1 rA'vvo 1 i" 1 i~\n y~\ iti't t 

\ukjL 1 clclLiUli cilt:! gy 


1 11 1 1 l-f 1 1 ^ /I oviTT T I /r\ r\\ Trio rrwyc^ 1 ti on oy m r 1 
l^Udi rVJ--'lVl cllcigy y/O Ul lllt: CUiicicxLlUll clitJl^y 1 


Be 


-14.556086 


-0.0527274 


-14.556123 (100.07) 


LiH 


-7.972557 


-0.0190867 


-7.9727078 (100.79) 


BH 


-25.058806 


-0.0569044 


-25.061771 (105.21) 


Li2 


-14.837571 


-0.0286889 


-14.839066 (105.21) 


BeH2 


-15.759498 


-0.0335151 


-15.761284 (105.33) 


H2O 


-75.735839 


-0.0546392 


-75.738582 (105.02) 


NH3 


-56.0586005 


-0.0693410 


-56.074805 (123.37) 
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Parameter |j. 



0.06 



0.05 



0.04 

to 

CD 

^ 0.03 

ro 

w 

a 

0.02 



0.01 



0^ 
-1 .446 



-1.445 -1.444 -1.443 -1.442 

Parameter ^ 



-1.441 -1.44 -1.439 
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FIG. 1: Cances et al., Journal of Chemical Physics. 



TABLE II: Cances et al., Journal of Chemical Physics. 



Oy o Lclll 




1 /^T"T"£i 1 Q T 1 /~»n on^ivmr 
^(Jl 1 CldLiUli CliClgjJ 


1 11 1 1 l-f 1 1 ^ /I oviTT T I /r\ r\\ Trio rrwyc^ 1 ti on oy m r 1 
l^Udi rVJ--'lVl cllcigy y/O Ul lllt: CUiicicxLlUll clitJl^y 1 


Be 


-14.613545 


-0.0467812 


-14.613653 (100.23) 


LiH 


-7.995678 


-0.0185565 


-7.9959693 (101.57) 


BH 


-25.171730 


-0.0630461 


-25.176736 (107.94) 


Li2 


-14.893607 


-0.0277581 


-14.895389 (106.42) 


BeHs 


-15.798440 


-0.0402691 


-15.801066 (106.52) 


H2O 


-76.120220 


-0.1401501 


-76.142125 (115.63) 


NH3 


-56.291315 


-0.1336141 


-56.318065 (120.02) 
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TABLE III: Cances et al., Journal of Chemical Physics. 



System Spatial basis size r CPU time (s) Newton iterations 

Be 5 25.7 2 

LiH 6 240.9 3 

H2O 7 958.8 4 

BeH2 7 1143.3 3 
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■107.4 



■107.6 



■107.8 



-108 



■108.2 



■108.4 



■108.6 



■108.8 




Hartree-Fock 
Full CI 

Dual RDM method 




1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 
N-N distance (Angstrom) 



